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1 Abstract 

Many genes have been identified as driving 
cellular differentiation, but because of their 
complex interactions, the understanding of 
their collective behaviour requires mathe- 
matical modelling. Intriguingly, it has been 
observed in numerous developmental con- 
texts, and particularly hematopoiesis, that 
genes regulating differentiation are initially 
co-expressed in progenitors despite their an- 
tagonism, before one is upregulated and 
others downregulated. 

We characterise conditions under which 
3 classes of generic "master regulatory net- 
works", modelled at the molecular level 
after experimentally-observed interactions 
(including bHLH protein dimerisation) , and 
including an arbitrary number of antago- 
nistic components, can behave " multi- 
switch", directing differentiation in an all- 



or-none fashion to a specific cell-type cho- 
sen among more than 2 possible outcomes. 
bHLH dimerisation networks can readily 
display coexistence of many antagonistic 
factors when competition is low (a sim- 
ple characterisation is derived). Decision- 
making can be forced by a transient in- 
crease in competition, which could cor- 
respond to some unexplained experimen- 
tal observations related to Id proteins; the 
speed of response varies with the initial con- 
ditions the network is subjected to, which 
could explain some aspects of cell behaviour 
upon reprogramming. 

The coexistence of antagonistic factors at 
low levels, early in the differentiation pro- 
cess or in pluripotent stem cells, could be 
an intrinsic property of the interaction be- 
tween those factors, not requiring a specific 
regulatory system. 

Abbreviations: bHLH, basic Helix-Loop- 
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Helix, Id, Inhibitor of Differentiation 

Keywords: multistationarity, cellular dif- 
ferentiation, cellular reprogramming, bHLH 
dimerization 



2 Introduction 

It has long been recognised that cellular 
differentiation could result from epigenetic 
memory, controlled by the dynamical 
properties of the same system, present 
with an i dentic al structure in all cells 
( DelbriickL Il949> ). rather than from a pro- 
gressive, irreversible loss of differentiation 
potential; a fundamental property of such 
a control system would be the presence 

of po sitive feedba c k cir c iiits Jxhomas , 

198lL jPlahte et all Hooj. ISnoussil. Il99^ . 



Gouzel Il998l , ICinauin and Demongeot 



2nn2[ ISou]eri2nml ). indeed, pioneer 
experiments showed that the genomes 
of some differentiated cell types retain 
the capa city to re genera te a whole or- 



ganis m 
192i), 



( Gurdon . 
and more 



1962L iGurdon et al 



recent experiments 



have strengthened the view that there is 
extensive plasticity in cell-fate determina- 
tion (r eviewed by Blau aiid Baltimorel . 



Blau and Blakelvl . ll999L and 



19911 

Theise and KrauseL 12002^ 



Bistable switches have been given 
a thorough th e oretic al investigation 



( Gherrv and Adleil . 1200(1 ). and have been 
constructed de novo flGardner et al 



2000[ l or modified ( Ozbudak et al 
There is evidence, discussed in 



section 

12.11 that cells undergoing differentiation 
sometimes face commitment decisions 
which involve more than two possible 
outcomes, but switches involving more 



than two variables have not been given 
extensive attention (we are not aware 
of any generic mathematical model that 
addresses cellular differentiation, with 
more than two possible outcomes). In 
the following, we discuss the relevance 
of these high-dimensional switches to 
the modeling of cellular differentiation, 
and investigate the properties of different 
molecular models, evaluating them with 
the current knowledge of the mechanisms 
of cellular differentiation. In particular, 
we test whether these models are able 
to display a coexistence of antagonistic 
factors at low levels, as decision-making 
with increased expression levels could be a 
relevant model of differentiation. 



2.1 Biological aspects 

2.1.1 Some commitments are irre- 
ducible to binary steps 

Cellular differentiation is often envisioned 
as a temporal cascade of decisions, by 
which cells restrict their potential fate fur- 
ther and further, until they reach a unique 
fate. It has been argued that each of 



these de cisions is binarv (|Brown et al 
19881 ISternberg and Horvit7j._ 11989 



Kaletta et al 



ml 



Lin et al 



19981 ). However, recent studies of 
hema topoeisis strongly sug gest other- 



wise fiRothenberg et al.l . Il999[ ). and point 



to models in which many cross-antagonising 
factors compete with each other (see be- 
low), receiving activation or inhibition 
from extracellular signals, leading to the 
progressive up-regulation of one specific 
factor, and down-regulation of all others. 
The hypothesis that decisions are more 
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complex than binary is also supported 
by the fact that the same cell type can 
be obtain ed by different develo pmental 



pathways (jRothenberg et al. 

Apart from hematopoiesis, two sys- 
tems have been described which seem to 
clearly involve a 3-outcome decision, irre- 
ducible to a sequence of 2 binary deci- 
sions: cells in the C. elegans hemaphrodite 
germline are directed to mitosis, differenti- 
ation as sperm, or differ entiation as oocyte 
( Ellis and Kimblel . lT995h . and founder cells 
of Drosophila mesoderm are directed to 
specific dorsal muscle or pericardial cell 
phenotypes by 3 m utually-repressive genes 



( Jada et all 1200 



Finally, in at least two instances of 
neural development, fate choices between 
a great diversity of possible outcomes 
have been shown, and are unlikely to 
be mediated by a series of binary com- 
mitments. This i s the case of olfac- 
tory development (ISerizawa et aL . 



Ebrahimi et al 



2000L 



which does not in- 



volve genetic rearrangements (jEggan et al, 



20041 ) , and of the regulation of hundreds of 



alternatively-spliced transcrip ts of a single 



gene in the Drosophila brain (jNeves et al, 



mi). 

Thus, it appears that model a, depicted 
in Figure is not the only possibility, and 
that model b of Figure [T] should also be 
taken into account. 

Having shown that high- dimensional 
switches are necessary for the mathematical 
modeling of some developmental decisions, 
we now turn to the way their structure 
should be modeled: differentiation factors 
are often antagonistic f section [2. 1.2|) . which 
doesn't prevent them from being sometimes 
coexpressed (section EUSI), and modulation 



a) Binary, hierarchic decisions model 
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Genes specific for cell-type 



b) Simultaneous decision model 




Figure 1: Arrows represent activation, 
and squares inhibition. Ad apted from 
Cinquin and Demongeot ( 20021 ). 
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of the interaction strength is a way differ- 
entiation is regulated f section [2.1 .41) . The 
basis for a mathematical formulation of the 
models is provided in section 1221 

2.1.2 Antagonism between differen- 
tiation factors 



Antagonism between genes driving differen- 
tiation to different fates has been repeatedly 
established; often, enforced expression of a 
differentiated phenotype, whether by spe- 
cific misexpression of a gene, or fusion of 
cells with different phenotypes, also leads 
to repression of the previous phenotype (re- 
pression of alternative fates has been pro- 
posed to be an essential mechanism of dif- 
ferentiation, reviewed by Corv . 19991 ). The 
idea of competition is reinforced by dose- 
dependency effects, shown for example by 
comparison of heterozygous and homozy- 
gous mutants, heterok aryon studi es, or 

knock-down mutation s ( Weintraubl . 1993 . 



repressive bHLH proteins (although not all 
possible cr oss-repressions have been char- 
acterised). iBriscoe et al.l ()2nnCll ) have also 
shown that a cross-repressive gene network 
reads out the Shh gradient in the neu- 
ral tube. Two sets of two cross-repressing 
genes have been identified, with a possi- 
bility that there is a larger, totally cross- 
repressive network (all the possible inter- 
actions do not seem to have been assessed 
yet). The competition can also happen 
by physical interaction between the fac- 
tors, rather than by cross-repression of 
transcription: in hematopoeisis, GATA- 
1, which drives eryt hroid and megakary- 
ocytic differentiation dKulessa et a . , 19951 
Visvader et al. . 19921 Iwasaki et al. . 2003[ ). 
and PU-1, a transcription factor essential 
for the expres sion of mvleo i d-spec ific genes 



(reviewed by Zhang et al. . 199fih. as we ll 



McDeyitt et al.l. I1997L review ed by lOrkinl . IZhang et al 



as B-cell specific genes (jChen et al. 
suppress each ot her's activity by phys- 
ical interac ti on (iRekhtman et ah , \l99i 



iNerlov et al 



Crittenden et all . 120021 1 , by monoal- This seems to be a general ph enomenon in 



lelic expres s ion o f a gene such as Pax5 
( Nutt et al. . I999I ). and by d osage effects of 



interacting bHLH proteins ()Zhuang et al 
Il99fi[ ). These effects argue that boolean 
models, in which a specific master gene 
would be turned on, initiate transcription 
of cell-type specific genes, and repress all 
other fates, are insufficient. 

Competition between cell-fate determin- 
ing factors has also been documented 
at the molecular level, for example in 
the case of neurogenesis, where bHLH 
proteins play a mai or role in specify- 



hematopoiesis (|Hu et 



19971 



1997L reviewed b 



Enver and Greavei 



I 



Cross and Enveil 
19981 ). 

In addition to repressing other genes, 
cell-fate determining factors often enhance 
their own expression; it has been proposed 
that this is a common property of "master 



switches" f Rothenberg et al. . 1999| ). 



ing neural subtype s dchien et ah j_ 19961 



Brunet and Ghysenl . |1 



Gowan et al, 



1)2001 ) have identified a network of 3 cross- 



2.1.3 Coexpression of antagonistic 
factors 

Coexpression of antagonistic genes has 
been shown both for closely-related 
lineages (for example coexpression of 
antagonistic hematopoeisis-related genes. 
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Miyamoto et a 



200 



2001 lAkashi et al 
^ reviewed by Orkin 



2003, lYe et al. 

200,4 transient prespore expression of a 
prestalk-specific gen e in Dictyostelium 



Termvn and Wiliamj . 1199.4 coexpression 
of lineage-spe cific genes in panc re as de- 



velopment, 
and coexp re ssion 



Chiang and Meltonl . 1200.1 




neurogenic genes 



Pierani et all |2mi| 



Briscoe et all I2000L although the latter 
may be due to a transient effect of the 
misexpression method), and between 
more distantly-related lineages (for ex- 
ample expression of neur al markers by 



Goolsbv et al. 



hematopoietic precursors 
200.1 . 

A semi- quantitative analysis of the ex- 
pression of ma ny hematopo ie tic ge nes was 
performed by Akashi et al.l ()2000f ). show- 
ing that lineage-specific (and antagonistic) 
genes were co-expressed at low levels in 
precursors, before respe ctive upregulatiori 
and downregulatio r i (see R.othenbergl . I2OO0I 
Zhu and Emerson . 200l for reviews). At 
an earlier stage of development, markers for 
different g erm lavers are also transien tlv co- 
expressed (jWardle and Smithl . l2004 ). 



2.1.4 Regulation of differentiation 

Some proteins have been shown to have reg- 
ulative effects on differentiation in many 
different cellular contexts, and would thus 
prove interesting to incorporate in models 
of cellular differentiation. 

• Id proteins, ubiquitously expressed 
during development, seem to act 
as inhibitors of cell differentia- 
tion, by sequestering ubiquitously- 
expressed class A bHLH proteins. 



preventing class B bHLH to form 



scriptionally 


active (Benezra et al.. 


1990l Garrell ann 


Modolelll. Il99n. 


Ellis et al., L 


1990[ 


reviewed by 


Norton et al.. 




I Norton. 2000l). 




Massari and Murre 
precise classification of HLH 
teins. Twist c an ac t in the same way 



( Spicer et al. . 1996t ). or in another. 



more d irect way, by binding to class 



MyoD (jHamamori et al.L Il997f l . 



Hesl, a bHLH protein, seems in 
many cases to be essential in the 
maint enance of an undiffer entiated 
state ( Kagevama et al. . I2OOOI ): its ef- 



fect can be mediated either by ac- 
tive repression, which involves the 
recruitment of Groucho, or by pas- 
sive repression, which involves hetero- 
dimerisation with other bHLH pro- 
teins. 

The PUF family of proteins re- 
presses the expression of many genes 
bv regulating t heir mRNA stability 



( Wickens et al.L 120021 ). and has been 
proposed to have the ancestral function 
of maintaining proliferation of stem 
cells; in C. elegans, sex-determination 
genes are regulated by PUF members. 



NF-kB has been shown to inhibit dif- 
ferentiation of mesenchymal cells, by 
destabilisation of the transcripts of 
Sox9 and MyoD, two transcription fac- 
tors involv ed in different di f ferent iation 
pathways (jSitcheran et al. [ 1200.1 . 
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All these differentiation-inhibiting pro- 
teins have a negative effect on the strength 
of transcription of genes which are essential 
in cell-fate determination. The models pre- 
sented below suggest that modulation of the 
transcription strength of proteins involved 
in cell-fate determination could allow for 
an initial co-existence of many antagonis- 
tic factors, followed by up-regulation of one 
factor at the expense of others, as the tran- 
scription strength is increased. 

2.2 Mathematical models 

The models studied here have an arbitrary 
number of components. Each variable rep- 
resents the intracellular concentration of 
a differentiation factor (called "switch el- 
ement" in the following), which enhances 
its own expression and represses that of 
all others (the system is symmetrical, in 
that any element has the same relation- 
ship with all others, and in that all ele- 
ments share a common set of parameters). 
The models can represent different forms 
of biological interactions. The terminology 
used below is that of transcriptional con- 
trol: each factor is supposed to be a pro- 
tein, which enhances the transcription of its 
own mRNA, and represses the transcription 
of the mRNAs for other switch elements, 
with or without physical interaction with 
other factors; clS db simplification, the trans- 
lation step is not taken into account in the 
model, and proteins are thus supposed to 
act directly on each other's concentrations. 
There is evidence that translational regu- 
l ation can play a major role in sorne case s 
( Wickens et 111 . \2QQ(i lOkano et all . 120021 ). 
In the following models, different forms of 
post-transcriptional control (by means of 



regulation of mRNA stability, or transla- 
tion of the proteins), can be represented 
in the same way as transcriptional con- 
trol. Downregulation of cytokine recep- 
tors has been observ ed prior to commitment 
( Kondo et alll2000h . and downregulation of 
receptors promoting expression of compet- 
ing factors could also be accounted for by 
the following models. 

3 kinds of models are studied below: 

• Mutual inhibition with autocatalysis: 
all switch elements repress one an- 
other, and enhance their own expres- 
sion. This is one of the simplest models 
one can think of that is able to achieve 
dominant expression of each of its ele- 
ments, depending on the initial condi- 
tions. 

• Mutual inhibition with autocatalysis, 
and leak: the same as the previous, 
with a supplementary term that repre- 
sents an identical, basal level of expres- 
sion, which is independent of any ele- 
ment of the network. This could corre- 
spond for example to a gene upstream 
in the differentiation hierarchy, which 
" primes" the lower level of the differen- 
tiation network, as has been proposed 
within the h ematopoietic d ifferentia- 
tion network jYe et all liooi reviewed 
bv lOrkinl . l20m 



• bHLH dimerisation: based on the 
class A / class B bHLH dimerisation dis- 
cussed above. 

The first two models can be viewed as 
a generic representation of the interactions 
between switch elements, while the third is 
based on an explicit assumption. All are 
formulated according to standard kinetics. 
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These models are cell-autonomous, and 
do not take into account "differentiation 
cues" that cells receive. The models could 
be extended to take into account either dif- 
ferent initial conditions, leading to various 
differentiated states, or different biases of 
the network (for example by providing a 
higher basal expression level of one of the 
factors). 

3 Results 

3.1 Mutual inhibition with 
autocatalysis 

Each switch element is supposed to undergo 
non-regulated degradation (modeled as ex- 
ponential decay, with an arbitrary speed 1), 
and transcription/translation with a rela- 
tive speed a. Each element positively auto- 
regulates itself, and represses expression of 
others, with a cooperativity c. Calling Xi 
the concentration of each switch element, 
the corresponding equations are, for n ele- 
ments: 



dXji 



-X\ + 



~l~ 



1 + S^ix^ 



GX, 



1 + S^ix^ 



The analysis is restricted to c > 1, as 
there is only one steady state (0) if c < 1. 
The results presented in appendix [X] show 
that it is possible for one switch element 
to be on and others off (for a > 2), but 
that if there is some cooperativity in the 



system [ie c > 1), it is impossible for more 
than 1 element to be on at the same time. 
On the contrary, if there is no cooperativity 
(c = 1), simulations show that a multitude 
of equilibria exist and are stable (switch el- 
ements in the "on" state can even coexist at 
different concentrations). Thus, the multi- 
stability behaviour of this system is tunable 
only by changes in the strength of the co- 
operativity, a mechanism which seems to be 
of unlikely biological relevance. 

3.2 Mutual inhibition with 
autocatalysis, and leak 

The model is the same as previously, except 
that each element has a " leaky" expression, 
modelled constant production term a. 
The equations become: 



GX\ 



-X\ 



1 + T.Ux\ 



1 + YIl^yLl 



+ a 



+ a 



The system is interesting only for c > 1 
(see appendix IB|i . If the leak is small, it 
doesn't have a major effect on the system, 
except when the cooperativity is close to 1: 
as shown in appendix El it is impossible for 
more than one switch element to be "on", 
at a much higher level than the leak level 
a. 

Even when the cooperativity is close to 1, 
it is still the case that only one variable at 
the same time can dominate all others; in 
order for that to happen, the transcription 
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strength must be sufficiently high. A simu- 
lation was performed for a cooperativity of 
1.1, with increasing a (see Figure |2)). All 
switch elements are initially coexpressed, 
and once a becomes sufficiently high, one 
switch element is upregulated, and others 
downregulated. 

The same pattern of coexpression fol- 
lowed by exclusive expression can be 
achieved with a decreasing leak (see Fig- 
ure O)), with the difference that the level of 
initial coexpression decreases slightly with 
time (this level is lower than the rela- 
tive maximum transcription strength a, but 
higher than the leak a). Once the leak has 
become sufficiently small, exclusive upregu- 
lation occurs. 

We show in appendix ^ that our models 
with mutual inhibition and autocatalysis, 
with or without leak, always converge to an 
equilibrium (and thus never oscillate). 

3.2.1 Effect of a perturbation 

If two switch elements are given initial val- 
ues close to the resting value one would have 
on its own, there is a transient drop in both 
values, until the higher one picks up and 
extinguishes the other (see Figure |3I). The 
initial drop is less pronounced than for the 
bHLH dimerisation model (see below). 

3.3 A model for bHLH pro- 
teins 

Each switch bHLH protein is supposed to 
bind to a common activator according to 
the law of mass action, with a binding con- 
stant K21 and a total quantity of activator 
a*. In turn, the heterodimer activates tran- 
scription of the corresponding switch pro- 




Figure 2: Time evolution of the concentra- 
tions of 4 switch elements (xi to X4), for 
the model with mutual inhibition with au- 
tocatalysis, and leak, with the transcription 
strength a being gradually increased over 
time. The 4 elements are initially coex- 
pressed at an identical level, which increases 
with a; when a reaches a threshold level, 
one element is upregulated, and others are 
downregulated. Parameters in the simula- 
tion were a = 2 and c = 1.1 Low, random 
noise was added to allow the system to es- 
cape the equilibrium as it became unstable. 
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1000 3000 3000 4000 5000 



Figure 3: Time evolution of the concentra- 
tions of 4 switch elements (xi to X4), for 
the model with mutual inhibition with au- 
tocatalysis, and leak, with the leak level a 
being gradually decreased over time. The 4 
elements are initially coexpressed at identi- 
cal levels (higher than the leak a because 
of autocatalysis); when the leak reaches a 
threshold level, one element is upregulated, 
and others are downregulated. Note that 
the scales for the Xi and for a are differ- 
ent by a factor of 11, equal to c/(c — 1) in 
this simulation. Thus, it is impossible for 
the curve of more than one Xi to be above 
that of a at equilibrium. Thus, in the boxed 
region, the system is in the process of re- 
sponding to the drop in a, and not at equi- 
librium. Parameters in the simulation were 
0" = 100 and c = 1.1 Low, random noise 
was added to allow the system to escape 
the equilibrium as it became unstable. 




2 4 6 8 10 



Figure 4: Time evolution of the concentra- 
tions of two switch elements (xi and X2), 
for the model with mutual inhibition with 
autocatalysis, and leak. The resting con- 
centration when one element is on and the 
other off is roughly 100. Initial concentra- 
tions differ by 10. Parameters are a = 100, 
a = 0.02, and c = 2. The trajectory is es- 
sentially the same for all a < 10, and very 
similar for initial concentrations differing by 
only 1. 
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tein only, with Hill kinetics and cooperativ- 
ity 2 (as with cooperativity 1, no interest- 
ing equilibria exist, as shown in appendix 
EJ- The equations are: 



dxi 



atxi 



= —Xi + a- 



dx 



n 



dt 



= -Xn + a 




These equations simplify to: 



dxi 
dt 



a 



^1 



aV^ + x: 



2 ' 

i 



with D = l + E^^^Xi, a, a = i^a/a? e K 
Parameter a is a measure of the harsh- 
ness of the competition between switch el- 
ements (it increases if K2 increases, ie if 
binding to the common activator occurs at 
higher concentrations, and if at diminishes, 
ie if the quantity of common activator goes 
down). The value of a is essential in de- 
termining the behaviour of the system. As 
shown in appendix O and summarised in 
section switch elements can coexist 

only if a is sufficiently low, ie if the com- 
petition is not too harsh (the lower a, the 
more switch elements can be non-0 at equi- 
librium). Figure El shows a simulation with 
a being increased over time; switch ele- 
ments are sharply turned off as a reaches 



threshold values. Figure IHl shows how an 
increase in a leaves only 1 switch element 
on, which remains exclusively on even if the 
competition level is relaxed to its original 
value. 

We prove in the appendix that the system 
always converges to an equilibrium for a > 
1/2; extensive simulations have also shown 



this to be the case for a < 1/2. 



3.3.1 Basins of attraction and times 
of response 

The cell fusion and reprogramming ex- 
periments, described below in section 14. 3| 
jvould lead to a situation where a switch el- 
ement, previously repressed, is brought to a 
level comparable to that of another switch 
element which was already expressed. This 
corresponds to an initial situation in which 
two elements are not at their resting value, 
which could also describe the situation in 
cells in the process of differentiating. For 
the models studied here, if 2 switch ele- 
ments are competing, 3 outcomes are pos- 
sible: the switch element at the higher con- 
centration completely represses the other, 
both coexist and reach a non-zero equilib- 
rium at the same value (only an element 
which started at the higher concentration 
can end up predominating), or both go to 
(extinction). Figures 171 to ITUl show which 
equilibrium the system converges to, clS db 
function of the initial state, for different val- 
ues of the competition parameter a (each 
domain from which the system converges to 
the same equilibrium is a "basin of attrac- 
tion"). If there are 3 switch elements com- 
peting, there are more possibilities, as 2 or 3 
elements can coexist in the "on" state. Fig- 
ures ^2 and El show the basins of attraction 
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Time 



Figure |): Same as Figureini but with a pulse 
of bife competition parameter a. The initial 
coji^jtions are such that the switch elements 
xist for low a; once a has sufficiently in- 
cr^cil^ed, only 1 switch element stays on, and 
dns on with all others off, even when a 
is'fitought back to its initial, low value. 



Figure 5: Time evolution of the concentra- 
tions of 4 switch elements [x\ to X4), in the 
bHLH dimerisation model, with competi- 
tion parameter a being gradually increased 
over time. The horizontal lines mark the 
values a = 1/4^, a = 1/3^, and a = 1/2^. 
Each time a reaches one of those thresh- 
old values, one of the switch elements is 
repressed. Low, random noise was added 
to allow the system to escape equilibria as 
they became unstable. In this simulation 
a = 100. 



of such a switch (the attraction basins be- 
long to the same system, but were split on 
two figures to prevent the outer ones from 
obscuring the inner ones). 

The speed at which the competition be- 
tween the switch elements is carried out 
could be crucial in determining the dynam- 
ical properties of differentiation. We thus 
computed the time it takes for the system 
to approach its equilibrium, starting from 
various initial concentrations of the switch 
elements (that time is colour-coded in Fig- 
ures m to IT^ . This time becomes dra- 
matically higher when the initial conditions 
come close to the borders of the basins of 
attraction (ze when concentrations are near 
a threshold around which the system con- 
verges to two or more different outcomes). 
The effect becomes even more pronounced 
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when 3, rather than 2, switch elements are 
competing (Figures ^2 and IT^ . 

To show the effect in more detail, in- 
dividual trajectories were plotted for a 2- 
dimensional switch fFiguresfT^andlT^ . For 
cell fusion and reprogramming experiments, 
the effect on the concentration of switch el- 
ements depends on the dynamics of nuclear 
import and export. Two types of initial 
conditions were used: two switch elements 
were given the concentration that one would 
have at rest, if it was "on" (Figure HH)), or 
two switch elements were given half that 
concentration (as cytoplasmic concentra- 
tions of proteins expressed exclusively in 1 
of 2 equally-sized cells should be cut by half 
upon fusion; Figure IT^. In both cases, the 
concentrations of the two switch elements, 
even for that which will eventually prevail, 
initially go down. The effect is more pro- 
nounced for higher values of the initial con- 
centrations, and for close initial values of 
the two concentrations. This could explain 
the transient extinction of expression of dif- 
ferentiated markers upon cell fusion (see 
Discussion): if there is sufficient coopera- 
tivity downstream of the switch element, its 
dip could be sufficient to provoke a tempo- 
rary extinction of expression of the proteins 
specific to the differentiated state. 



3.3.2 Extinction domain 

2 

For a > Tr^-TTT, there is an extinction do- 
main around the diagonal Xi = .. = 
Simulations show that the domain is very 
restricted until a becomes very close to its 
upper threshold value, at which non-0 equi- 
libria cease to exist (see Figures IHl and imjl . 




Figure 7: Colour-coded time of convergence 
(as defined in Appendix ID. 2j) . function 
of the initial conditions in xi and X2. From 
the initial conditions to the left of the red 
region, the system converges to X2 on and 
Xi off, and the opposite for the initial con- 
ditions to the right of the red region. Pa- 
rameters are a = 0.4 and a = 100. xi and 
X2 range from to 300. 
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Figure 8: Same as Figure [3 but with a 
lower value of a, giving a large domain of 
convergence to an equilibrium where Xi and 
X2 coexist. Domains of convergence are in- 
dicated, and are separated by the two yel- 
low lines. Parameters are a = 0.1 and 
a = 100. Xi and X2 range from to 300. 




Figure 9: Same as Figure [3 but with a 
markedly higher value of a. There are two 
main domains of convergence (to one switch 
element on and the other off), and a third 
domain of convergence to (complete ex- 
tinction of the switch), in a region very close 
to the upper part of the diagonal (for clar- 
ity reasons, the region is indicated as larger 
than it actually is). Parameters are a = 15 
and cr = 100. Xi and X2 range from to 
300. 
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Figure 10: Same as Figure [71 with a close 
to the threshold above which is the only 
equilibrium. The region from which the sys- 
tem converges to has expanded. Parame- 
ters are a = 24.75 and a = 100. xi and X2 
range from to 300. 

3.3.3 Summary of a threshold values 

For the system to be able to sustain k switch 
elements "on" at the same time, the condi- 
tion a < l/k"^ must be met (for a ^ 1, 
this condition is also sufficient). Thus, for 
a > 1/4, only 1 switch element can be on 
at a time. The corresponding equilibrium 
value is a decreasing function of a. For 



« > Tr^-TTT, there is an "extinction do- 
main" around the diagonal xi = .. = Xn'- 
matching sufficiently closely the concentra- 
tions of the switch elements, at whatever 
value, makes the system switch off all switch 
elements. For large a, the extent of this do- 
main is small, except in a very narrow range 



of a values. Finally, for a > -^^^z^, a con- 
dition which becomes a > a/4 for large a, 
there are no non-0 steady states. 




[ 

■-,03 



-^■ll = i], iji = (I 



Figure 11: Times of convergence as a 
function of the initial condition, for a 
3-dimensional switch. 4 unconnected 
domains of convergence to the same 
equilibrium are shown. For visibility, 
the 3 other domains are shown in Fig- 
ure ^1 Parameters are a = 0.1 and 
a = 25. A rotation movie is available at 



http: / / www-timc.imag.fr/01ivier.Cinquin/High-dimensionaL 
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Figure 12: Domains in which the 
same switch as in Figure ^2 converges 
to a state with 2 switch elements 
A rotation movie is available at 



Figure 13: Time evolution of the concentra- 
http://www-timc.imagir/Qlivier.Cinquin/Higfc4ip^i|Lyi3na.l^^i^yfe 



on. 



concentration when one element is on and 
the other off is roughly 8. Initial concentra- 
tions differ by 0.7 (a), or 0.1 (b). Notice the 
differences in scales. Parameters are a = 50 
and a = 500. 
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4 Discussion 

4.1 Co-expression properties 

Of the models proposed here, if the cooper- 
ativity of activation is considered to be con- 
stant, only the model with bHLH dimeri- 
sation is capable both of exclusive expres- 
sion of an arbitrary number of switch el- 
ements, and coexpression at similar levels 
of all elements. This behaviour is finely 
tunable with the key competition param- 
eter deriving from the availability of the 
bHLH hetero-dimerisation partner, with 
lower availability being unfavourable to co- 
5 10 1 5 20 expression of the antagonistic factors (see 

below for a further discussion). 

The model with mutual inhibition, auto- 
catalysis, and leak can express no more than 
one switch element at a level higher than 
the other ones, and is thus less amenable 
to progressive elimination of unwanted fac- 
tors in the course of differentiation. In order 
for coexpression to occur at a significantly- 
higher level than the leak, the cooperativ- 
ity of the system must be close to 1. If 
differentiation was controlled by a network 
of this kind, initial coexpression could be 
maintained either by a low transcriptional 
strength in the system (which is consistent 
Figure 14: Same as Figure CHI but with ini- with antagonistic factors being expressed at 
tial concentrations at roughly half the equi- a lower level in the un-differentiated state), 
hbrium value when one element is on and or, as has been suggested, by regulated 
all others off. Initial concentrations differ degradation of mRNAs. 
by 0.5 (a), or 0.1 (b). Interestingly, the multistability be- 

haviours of a switch based on bHLH-like 
dimerisation and that of a switch based 
on direct cross-repression are qualitatively 
different: the former can maintain many 
variables on at an equilibrium only if those 
variables are sufficiently high (compared 
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to the transcription strength), while the 
reverse is true of the latter. 

We previously studied networks of cross- 
repressing factors, in which the factors 
do not enhance their own expression 
( Cinquin and Demongeot . I2nn2t ). We did 



not include this kind of model in the present 
study, because for one factor to be able to 
dominate all the others, it had to be as- 
sumed that the cooperativity of the network 
was very high, an assumption which is pos- 
sibly not realistic. 



4.2 Peaks of differentiation in- 
hibitors 

According to the paradigm of inhibition 
of differentiation by sequestration of class 
A bHLH proteins, the quantity of binding 
partner should be low prior to differentia- 
tion, and the competition parameter a in- 
troduced earlier should thus be high. Re- 
lieving inhibition of differentiation, by in- 
creasing the quantity of binding partner, 
and thus decreasing competition, cannot 
move the bHLH dimerisation network from 
a state where it supports coexpression of 
many switch elements, to a state where only 
one is expressed. Also, increasing the tran- 
scription strength of the network a does not 
destabilise equilibria. 

It is thus impossible to account for the 
selection of one differentiation outcome 
by increasing the availability of class A 
proteins (for example by downregulation of 
Id proteins). However, it is still possible 
that the competition strength, even in the 
undifferentiated state, is sufficiently low for 
many switch elements to be co-expressed. 
A potential mechanism to select 1 element 



and extinguish all others is then to tran- 
siently increase the competition strength, 
for example by transient up-regulation of 
Id proteins, just prior to differentiation 
commitment (a corresponding simulation 
is shown in Figure IH))- This is in fact 
what has been experimentally observed 
on independent occasions, on a short time 
scale, durin g in vitro d ifferen tiation of 

E 



osteoblasts (iQgata et al 



( Nagata and Todokoro 



( Andres-Barauin et al.l. Il997^ 



cells fIStewart et a) 
(ILanglands et al 



cells fllshiguro et~ . 1996[). 



I993II , neurons 
19M), myeloid 



flHouldsworth et al 



Chambers et al 



Em 



2D 



astrocytes 
Schwann 



keratinocytes 
germ cells 
and fibroblasts 



_ 2nnih 

20031 ). No rationale for 
these transient effects had been proposed 
so far. 

When Id proteins are not up-regulated, 
other proteins could play the same role of 
increasing competition in the bHLH net- 
work. For example, Hes-1, which also has 



class A sequestering activity (jSasai et al 
1992[ l. is transiently upregulated upon dif- 
ferentiation of an immortalised hair cell 



line (|Rivolta et al.l.l2002h. an immortalised 



Ohtsuka et all Il998h 



neural cell line 
neuroblastoma (iGrynfeld et al.L |2| 



and 
(al- 
though its role in the latter case is compli- 
cated by the fact that it also binds Id pro- 
teins and i nterferes w' i th Id2 /E2-2 complex 



formation, iJogi et al.l . l2002h : the transient 



Hes-1 expression is concomitant with down- 
regulation of the bHLH protein HASH-1. 
Hes genes are dominant represso r s wit h 
many targets ( Barolo and Leviii3 . 1997), 
and could also directly repress many ele- 
ments of the network, which can be modeled 
by a decrease in the transcription strength 
0", and has the same effect of destabilising 
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equilibria where many elements are coex- 
pressed. 

Finally, erythroid-specific genes have 
been observed to be transiently downreg- 
ulate d upon induced, i n-vitro differentia- 



tion (jLister et all Il995( ). which could also 
be explained by transiently-increased com- 
petition in a bHLH dimerisation network. 

4.3 Dynamical properties 

Analysis of the proposed dynamical sys- 
tems shows that the time to convergence 
can widely depend on the initial condition. 
Convergence can be relatively very slow 
when initial conditions are near a threshold 
around which the system converges to two 
or more different outcomes. This is for ex- 
ample the case when 2 or more " switch ele- 
ments" are at roughly equal concentrations, 
higher than that of others (the more ele- 
ments in competition, the slower the com- 
petition becomes). It is interesting to note 
that slow effects are observed in induced- 
transdifferentiation experiments, and in cell 
fusion experiments. 

• Fibroblasts reprogrammed to T-cell- 
like cells need to be incubated for many 
days before they acq uire detectable 



T-cell characteristics fflakelien et al 



2002| ). This may be due to the fact that 



fibroblast master genes are expressed 
at a high level, and the counter-acting 
T-cell master genes, introduced by per- 
meabilisation of the membranes, are 
also present at a high concentration. 
An effect of the relative levels of cy- 
toplasmic factor concentrations could 
be tested by incubation in T-cell and 
fibroblast cytoplasmic extracts, mixed 



at different ratios. Further investiga- 
tion of master networks could involve 
incubation of cells in cytoplasmic ex- 
tracts of 3 or more cells-types (or tran- 
sient misexpression, at controlled lev- 
els, of antagonistic master genes). 

In hepatoma-fibroblasts hybrids, ex- 
tinction of al bumin production can 
take days ( Mevel-Ninio and Weissl . 
1981[ ). Most interestingly, some hy- 
brids show reexpression of albumin 
after extinction. These two outcomes 
can be accounted for by the models 
proposed above: when two antagonis- 
tic "switch elements" are coexpressed 
at a high level (which probably cor- 
responds to the fusion experiments, 
as upon fusion the protein contents 
of the cells, which are of different 
phenotypes, are mixed), it is possible 
for the system to revert to a state 
where all switch elements are turned 
off (total extinction), or for the two 
switch elements to decrease to a low 
level, before the trajectory of one 
of them picks up and goes back to 
a high state (extinction followed by 
reexpression) . 

Activation of the myogenic phenotype 
also takes place on the scale of days, 
when muscle cells are fused to var- 
ious other cell types, a delay which 
was suggested not to be lin k ed to 



DN A duplication jBlau et al. 
Blau and Blakelvl . ll999L for 



see 



tensive review) 



1985[ 



an ex- 



Also, it could be that the progres- 
sive upregulation of differentiation-related 
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genes observed during hematopoietic devel- 
opment is a cell-autonomous consequence of 
the slow dynamics of a switch network. 

4.4 Stochastic outcomes 

It has been observed in the studies cited 
above that heterokaryons with the same 
number of nuclei coming from each donor 
can have different differe ntiation responses. 
Blau and Blakelv ( IQQfli l suggested that the 
differentiation outcomes are regulated by 
a network of opposing factors. Within 
this framework, stochastic responses to cell 
fusion can readily be explained by slight 
differences in the quantities of differentia- 
tion factors contributed by each cell type, 
which tip the balance one way or the other. 
The network determining cell fate would be 
most sensitive to random noise when the 
factors are in roughly equal concentrations. 
The sensitivity to molecular noise of the 
networks proposed here would be interest- 
ing to study, as it has been proposed that 
some types of differentiation during devel- 
opment could have a stoch astic aspect ffor 



exam ple in adipoEfene sis. iTchkonia et al 



2002 , or hematopoiesis. lEnver and Greaves 

mm- 



4.5 



Evolvability of switch net- 
works 



In addition to having a coexpression be- 
haviour which is easily tunable by one pa- 
rameter, the generic bHLH network stud- 
ied here has the advantage of being per- 
haps more easily evolvable than a network 
in which every element actively represses 
all others: the interaction needs only take 
place between every element and a common 



activator (which requires n interactions, in- 
stead of n (n — 1) /2). bHLH networks have 
been suggested to evolve mainlv bv single 



gene duplication events (jAmoutzias et al 



20041 1 ■ maintaining a topology in which ev- 
ery element of the network interacts with a 
restricted number of "hubs". 



5 Conclusion 

The models presented here could be useful 
in understanding development, as well as 
cell-fate reprogramming (which can be in- 
duced artificially, b ut has also been shown 
to happen 

We have derived general results about the 
dynamics and co-expression properties of 
switch networks, and shown the flexibility 
of bHLH dimerisation networks. Of the 
networks studied here, which were gener- 
ically formulated with usual kinetic equa- 
tions, only a subset can co-express antag- 
onistic elements at a similar level, higher 
than the basal level: those with mutual 
inhibition, autocatalysis, and leak (but 
only when the cooperativity is very close 
to 1, and the transcription strength suffi- 
ciently low), and bHLH dimerisation net- 
works (when the competition is sufficiently 
weak). This restricts the classes of mod- 
els which can reproduce experimentally- 
observed co-expression of antagonistic fac- 
tors, as well as showing how it can occur. 

Strikingly, even though bHLH networks 
are the most apt to coexpression of antago- 
nistic elements, the selection of one element 
requires a transient increase in competition, 
which is not what is thought to happen over 
a long time scale in the course of differen- 
tiation. Transient, hitherto- unexplained in- 
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creases in competition have however been 
observed in a few cell lines upon differenti- 
ation, and could be a general phenomenon. 

In order to model specific differentiation 
events, these networks would probably need 
to be extended to take into account combi- 
natorial interactions, which could compli- 
cate their behaviour. The models would 
also gain from being extended to take 
into account non-symmetrical networks, in 
which some switch elements are stronger 
than others, and stochastic kinetics. 
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an equilibrium, from any initial condition, 
will be derived in a more general context, 
in section IB. II In the rest of the appendix 
we assume c > 1. 

A. 2 One on, all others off 
A. 2.1 Equilibrium existence 

The steady-state equations are 

Vj,%(l + SILiX^)=a4 



le 



xf ^ = ; (1 + SILiX?) or X, = 



Re-arranging the first equation. 



^1 + Yji^jX 



Let /(x) = -x'^ — . Then fix) = 
fx'^-i - (c - l)V-2. f (x) < iff fx < 
c — 1. The minimal value of / over the 
positive real set is /(—a) = - [—aY — 

The equilibria studied here are such that 
only Xj is non-0, for some j. There are ei- 
ther or 2 solutions, 2 iff 
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c-l 



- > - 

c a 



a > c 



c-l 



Incr > 



c-l 
lnc + ln(^) 



c-l 



(4) 



In {-^Y " is an increasing function of c, 

and limc^ooln (^j"^ ^ = 1- — is decreas- 
ing for c > e ~ 2.7. The right-hand side 



of equation |3] has a maximum for c = 2, 
of about 0.7, matched by a = 2. Thus, for 
a > 2, there are two equilibria. Both large c 
and large a are favourable to the existence 
of an equilibrium with one variable domi- 
nating all others. 

A. 2. 2 Local stability analysis 

It is useful, for the Jacobian term computa- 
tions to follow in the rest of the appendix, to 
note that if ^(x) = g'{x) = f^^- 

If Xj is at a non-zero steady-state and 
Vz 7^ j, Xi = 0, and if c > 1, the stability 
at that steady state depends only on the 
sign of the (j, j) coefficient of the Jacobian 
matrix (this coefficient will be called Jjj in 
the remainder of the appendix). 
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with 
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at equilibrium 
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axj 

the equilibrium is stable iff 



> -, ie 1 x^ > c 



a 



It is possible to give a sufficient condi- 
tion for the equilibrium with the greatest 
solution to equation IHl to be stable. Let 
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fix) = X'- ax'-\ If / ((f) 

then the greatest root of equation [HI will be 

1 

greater than (|:) "'^ , and the corresponding 
equilibrium will be stable. A sufficient sta- 
bility condition is thus 



< c- 1 



Numerical investigation shows that this 
condition is met for cr > 2. 

A. 3 k variables on, others off 

With identical parameters, there can be no 
equilibrium with 2 variables having differ- 
ent, non-zero values. 

At any equilibrium, variables can be 
renumbered so that, in the Jacobian ma- 
trix, variables at form an independent 
block. This block is stable, and the sta- 
bility of the whole system depends only on 
the block formed by non-0 variables. Thus, 
in the following we suppose that no steady- 
state variable has for a value. 

For i ^ j. 
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With the same kind of ana lysis as in 
Cinquin and Demongeot ( 2002t ). the equi- 
librium is stable only if 
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With the definition of the equilibrium. 
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(Tc{l + {k- 1) 



X ) X 



c-1 



x'K ^x'^^ - {1 + (k - 1) x") 



a 



-x"'^ > kx" + 1 
c 

Again with the definition of the equilib- 
rium, 

-x"-' > x'-\ 
c 

ie c < 1, in which case no interesting equi- 
libria exist. 



B Analysis of mutual in- 
hibition with auto- 
catalysis, and leak 



If a > 0, c > 1, and one of these inequalities 
is strict, the function /(x) = x^~'^ — ax~'^ 
can take the same value for at most 2 pos- 
itive values of x. Thus, there are only two 
values a variable can take at a given steady 
state (0 cannot be a steady state value). If 
two different equilibrium values are taken 
by some variables, one of these values is 
higher than a-^-, and the other lower. 

If a > and c = 1, the system only has 
one equilibrium, with all variables equal. 

B.l Convergence 

Let Ui = ^/xi, and 
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Xi 



2m 



2c- 1 



1 + SLiZ/. 



2c 



+ - 



Thus, P is a potential for the system. 

If its equihbria are isolated, a gradient 
system converges to a steady-state regard- 
less of the initial conditions. It is shown 
below that the number of solutions of the 
system is finite when the cooperativity c is 
an integer, and the system thus always con- 
verges to a steady state (we expect this re- 
sult to also hold for non-integer values of 
c). The model without leak corresponds to 
a = 0, and this convergence result thus also 
applies to it, for c > 1. 



nx 



c+l 



(a + na) x'^ + x, 



If fix) 

f'{x) = (c + l)nx'' — c((T + nQ.)x'^~^ + 1, 
f"{x) = c{c+l)nx''~^-c{c-l){(r+na)x''~^. 

f" f;r77XTT + '^")') = 0- /' takes negative 



^rt(c+l) 

values iff /' 



n(c+l) 



a + na) I < 0, which is 



a necessary condition for the existence of 3 
equilibria with all variables on. 

The dynamics of the system constrained 
to Vz, defined by 



ax 



X 



-X 



a 



1 -|- nx'^ 

The sign of x{t) is the opposite of that of 
f{x(t)). Because x„ is such that /'(xj) < 0, 
it is easy to see that the steady state x„ 
is unstable for the constrained system, and 
thus for the full system. 

B.2.2 Local stability analysis 

With a leak a, equation IA.3I becomes 



x" < 



a X 



c+l 



cix- a) 



B.2 Steady-state analysis: all 
at the same value 

B.2.1 Equilibrium existence 

V J, {x, - a) (1 + T^UxX) = ax'^^ 
If V j, Xj = X, 



nx 



x'^^^ — (a + na) x'^ + x — a = 



There is at least one solution, maybe 3 
(or 2 in degenerate cases) depending on the 
parameters. The solutions are noted xi, Xu, 
and Xh, with xi < Xi < Xh- 



l + kx" < 



c (x - a) 



< 



a x"+i 



X — a 



c [x 



a] 



X — a < -X 
c 

Thus the stability condition lA.3l is met iff 
X < Ci-;rr[ (in thai case, since non-diagonal 
terms of the Jacobian are obviously nega- 
tive, diagonal terms are also negative, and 
the equilibrium is stable). Since solutions to 
equation |H1 can be made arbitrarily high by 
increasing a, increasing a past a threshold 
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value (other parameters being equal) will 
prevent the existence of a stable equilibrium 
with all variables equal. 

B.3 k on, k < n 

Let p = n — k. 



{xi — a) (1 + px'^ + kxl) = ax 



px'i^^ 



{pa + cr) 



kxt' 



a (1 + kx 
— {ka + 0") 
ail+px'i) 



(1 + kxl) xi 


(1 +px'i) Xh 





variables having value Xi is strictly greater 
than 1. 

Thus there are only two possible kinds 
of stable equilibria: all variables equal, in 
which case the equilibrium value is lower 
than a-S", or one higher than all the other 

c— 1 ' ° 

ones (in which case the lower ones are lower 
than, and the higher one greather than 



a 



c-\r 



C Analysis of the bHLH 
model 

Without cooperativity in transcriptional 
activation by the bHLH dimer, there is only 
one stable steady-state: 



Choosing for example the graded lexico- 
grap hic order o v er Cb /,^^], theorem 5.3.6 
from ICox et al.l ( 19961 ) shows that the sys- 
tem has a finite number of solutions, when 
c is an integer. 

We have 



Ji.i = -1 + caxl 



J. . ^ _r^^c-i^ 



If 



Ji i Ji 



c-1 ^ 



-1 + c- 



Consider the reordered Jacobian matrix, 
with k variables "on" with a value Xh-, and 
p "off" with a value xi {k + p = n). 

It follows from the analysis in section lU^I 
that the equilibrium can be stable only if 
Ji,i—Ji,j < (ie Xi < a^), if the number of 



Xi Xi 



■1 + 



cr 



a 



(1 + S^=,x,) 



Xi 



If at some steady state k variables are on 
and share a common value x (variables at 
a steady state, if not 0, share a common 
value) , 



kax + X + a 



X 



a — a 



ka + 1' 



and if Xp(to) = 0, 



a 



J, 



kax + a 
a — a 



a {ka + 1) 



>0, 



and Jp^i = for p ^ I, proving the unsta- 
bility of the steady state. 
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In the following, it is assumed that tran- 
scriptional activation occurs with coopera- 
tivity 2, and the steady-state equations be- 



come 



X7 



Xi 



axi 



(9) 



C.l Dynamical analysis 



is a stable steady state. If Xj(0) = 0, 
then V t > 0, Xi{t) = 0. If Xi{0) > 0, then 
V t > 0, Xi(t) > 0. One can thus suppose 
that \/ i,\f t > 0, Xi{t) > 0. Consider a 
state in which there is one variable strictly 
superior to all others (ie, a state not belong- 
ing to the line Xi = X2 = ■■ = Xn)- Suppose 
without loss of generality that the variable 
in question is xi. Consider the function 



/i(^) 



x; 



aD"^ + xi 



fi{x) = 2aDxi 



X- 



{D - xi) - xjll 



X 



2^2 



Suppose that (t/i(0) > x(0). In this case, 
/i (0) > 0, and Xi and / are strictly increas- 
ing functions of time. If o"/i(0) < a;i(0), 
then /i (0) can be negative or positive. In 
the first case, xi is decreasing as long as /i 
is. If at some time to ^fiiio) ^ ^(^o); then 
for t > tQ, Xi and /i are increasing func- 
tions of time. Thus there can be at most 
one change in the monotony of Xi. Thus 
limt^oo exists. Since x'l exists and is 
bounded on any trajectory, limt^oo ^Ji (t) = 
0. All trajectories thus converge to a steady 
state where Vj > 1, Xj = Xi or Xj = 0. 

If 3 t st V 72 > j > 1, xi{t) = Xj{t), the 
system is brought back to one dimension. 
Note that it is impossible for any variable 
to outgrow Xi. 

C.2 Steady- state analysis: 
variables on at the same 
value 

C.2.1 Equilibrium existence 

Variables zero at the steady state can be 
discarded from the analysis. If k variables 
are non-0, and are all equal, to x 7^ 0, 



(aD^ + xlf 



fi{x) = xr 



2aDxi 

xiXi (xi — Xi) {aD"^ — xiXi] 



{aD^ + xj) {aD^ + , 

For a > 1/2, the second term is positive. 
We have 

dxi(t) 

= (^fi{x) -xi 

We first consider the case in which V t > 
0, V n > j > 1, xi > Xj. 



x^{l + Pa)+x{2ka-a)+a = (10) 
Solutions are 



a 



- 2ka ± y/a"^ - 4a (1 + ka) 
2{l + k^a) 

A sufficient and necessary condition for 
the existence is 

ka + 1 



4q!- 



< 1 
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It will be shown below that, at a stable 
steady- state, there is at most 1 non-0 vari- 
able which can be different from other non-0 
variables. If there is such a variable, equal 
to the equation for the value of other vari- 
ables becomes 



Using the steady state equation IHl 

Ji,i = -l + — {D{D~ X,)) = -1 + 
axi 

2 

— {a — Xi — aD) 

(7 



x^ (l + k^a) + X {2ka {I + y) - a) + 



+ = (11) 



Solutions are 



Jii = \ {xi + aD) = 1 (a -|- Xj (1 -|- ka)) 

a a 

The diagonal terms are negative for 
a/2 - a 



o - 2ka {l + y)± a'^ - 4a (1 + A;a y) (1 + y) 



Xi > 



ka 



2(1 + Fa) 
and the condition for a solution to exist 



"The off-diagonal terms are given by 
2 2axj + 2a {D 



4a (1 



ka + l + y 

y) — < 1 



-ax. 



a" 



The solutions for y are 



J, 



-2crax, ■ 



(aD2 + 
D 



Ji 



{aD^ 
2a 



X 



2^2 



D 



a 



a - 2a (1 + kx) ± ^Ja'^ - 4a (1 + a + kx) (1 + kx) 
2(1 + a) 

C.2.2 Local stability analysis 

Thus, a necessary condition for the equi- 
Variables zero at the steady state can be ^-^^-^^ ^^^^^^ 

discarded from the analysis. 

Using 



-1 + 2^ 
a 



V Xi st Xi 7^ 0, Xi> cr/2 



(12) 



fex^ + 2cx 
ax2 + hx + c (ax2 + hx + c 



2 ' 



one derives the diagonal term of the Ja- 
cobian (with h = 2a{D — Xi) and c = 
a{D — Xj)^): 



Ji. 



— 1 + 2a ax. 



D{D 



Xi\ 



\aD^' 



x; 



2^2 



This is possible if and only if a < l//c^ 
and a > 2^^. 

Condition^Jis stronger than the require- 
ment for the diagonal element to be nega- 
tive (and is thus also a sufficient condition), 
and can never be met by variables equal to 
the lower solution of equations ^1 or • 

Thus, for any value of the transcription 
strength a and for any number of coexis- 
tant variables k, sufficiently low values of 
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ot make the equilibrium stable. If there is a 
stable equilibrium with k variables on, there 
is also a stable equilibrium with p variables 
on, for 1 < p < fc. For sufficiently large a, 
the necessary condition a < becomes 
sufficient for stability (see FigureElfor an il- 
lustration of the validity of this condition). 

C.3 On at different values 

If at steady state, Xi 7^ Xj and both are 
non-0, then 



condition for stability of an equilibrium is 
e > b and c> a. In particular, there can be 
at most 1 variable having value. 



x1 — axi = x| — axj (= —aD^) 

There are thus only two possible non-0 
steady-state values, noted Xa and Xb, with 
Xa < Xb- Noting P{x) = x^ — ax, and sup- 
posing that Xa and Xb exist, P'{xa) < 0, ie 

cr 

Consider the Jacobian matrix of the sys- 
tem, reordered so that variables having Xa 
as a value come before those having 
value: 



k 



c ■■■ c a fi /i 



V/2 /2 e ■■■ e bj 

With the appropriate eigenvectors, it is 
easy to show that b — e and a — c are eigen- 
values for this matrix, of order k — 1 and 
p — 1. Thus, if /c > 1 and p > 1, a. necessary 
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More precisely, the characteristic polynomial of the matrix is 

P{x) = (a — c — x)''^^ (6 — e — xY~^ [x'^ — x {a + b + {k — 1) c + {p — 1) e) 
+ {p-l)ea+{k- 1) cb+{k- 1) {p -l)ec + ab- kpf 1/2] 

Suppose thus that the number of variables having values Xa is 1. Then, a sufficient 
condition for instability of the equilibrium is 

{p -l)ea + ab- p/1/2 < 

Notice that in this case /i = /2 = e. The sufficient condition for instability can thus be 
written 

e (pe — {p — 1) a) — ab > 
Replacing with the equilibrium values. 



-2aD f -2aD , / 2 , 

P^_ {p-1) (1- -{xa + aD) 



a 



a 



a 



l-^ixa + aD)] ( 1 - ^ {xb + aDU > 



l--{xa + aD)j Up- 1) ^— -l + -{xb + aD)j + pi —— \ > 



^ 2 , f 2aD ^ 2xb 
l--{xa + aD)] p 1 + — 



(2aD\ ^ 
p \ I > 



■) 



2aD ( 2xa\ (2xb \ ( 2 , 
P (1 ^ 1 + ( — - 1 ) ( 1 - - (x, + a/}) ) > 



a 



a 



a 



2Xn 



cr 



2aD 2xb 

p 1 1 

cr a 



a 

2aD (2xb 



- 



2aD ( ^ 2 , 

P + 1 ypXa + Xb, 

o \ a 



a \ cr 

23^a \ ( 2Xb 



a 



1 > 



1 > 



a 
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The first term is positive because the val- 
ues of Xa and Xi, are symmetrical with re- 
spect to cr/2. The second term is also posi- 
tive, and the sufficient condition for the in- 
stability of the equilibrium is thus met. 

Thus, there is no stable equilibrium with 
non-0 variables having different values. 



D Methods 

D.l Numerical integration 



results). The stepsize of the Runge-Kutta 
algorithm was kept lower than 0.3. 

D.3 Simulations with time- 
dependent parameters 

In order for the system to leave steady 
states which had become unstable because 
of changed parameters, small random per- 
turbations were applied (each variable was 
multiplied by a random number uniformly 
chosen in [0.99 .. 1.01] every 30 time units). 



All integration was performed with a 
custom-written implementation of the 
4th-order a daptat i ye ste psize Runge-Kutta 
algorithm ( Pressl . Il992[ ). with 10~^ rela- 
tive accuracy. Source code is available at 

|http : //www-time . imag.fr/01ivier . Cinquin/ada/ ada_blas_runge_kutta . html 

The data was plotted using GMV or gnu- 
plot. 



D.2 Computation of conver- 
gence times 

A custom program was written to do the 
following, starting from a regular 200*200 
grid of initial conditions (for 2D systems), 
or a 50*50*50 grid (for 3D systems), with 
Vi 7^ j, Xi 7^ Xj, to avoid reaching unsta- 
ble steady-states: (1) integrate the system 
until a steady-state is reached (as defined 
by the sum of the absolute values of the 
derivative vector elements begin lower than 
10~^) (2) start the integration again, with 
the same initial conditions, stopping when 
the system gets close enough to the previous 
steady-state (each variable with 10% of its 
steady-state value if it's not 0, lower than 
0.15 if it is 0; moderate changes in these ar- 
bitrary values do not significantly affect the 
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